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One of the goals of this program was to develop the experimental and analytical/computational tools 
required to predict the flow of non-Newtonian fluids through the various system components of a 
propulsion system: pipes, valves, pumps etc. To achieve this goal we selected to augment the capabilities 
of NASA’s Generalized Fluid System Simulation Program (GFSSP) software. GFSSP is a general- 
purpose computer program designed to calculate steady state and transient pressure and flow distributions 
in a complex fluid network. While the current version of the GFSSP code is able to handle various 
systems components the implicit assumption in the code is that the fluids in the system are Newtonian. 
To extend the capability of the code to non-Newtonian fluids, such as silica gelled fuels and oxidizers, 
modifications to the momentum equations of the code have been performed. We have successfully 
implemented in GFSSP flow equations for fluids with power law behavior. 

The implementation of the power law fluid behavior into the GFSSP code depends on knowledge of the 
two fluid coefficients, n and K. The determination of these parameters for the silica gels used in this 
program was performed experimentally. The n and K parameters for silica water gels were determined 
experimentally at CFDRC’s Special Projects Laboratory, with a constant shear rate capillary viscometer. 

Batches of 8:1 (by weight) water-silica gel were mixed using CFDRC’s 10-gallon gelled propellant 
mixer. Prior to testing the gel was allowed to rest in the rheometer tank for at least twelve hours to ensure 
that the delicate structure of the gel had sufficient time to reform. During the tests silica gel was pressure 
fed and discharged through stainless steel pipes ranging from 1” to 36” in length and three diameters; 
0.023”, 0.032”, and 0.047”. The data collected in these tests included pressure at tube entrance and 
volumetric flowrate. From these data the uncorrected shear rate, shear stress, residence time, and 
viscosity were evaluated using formulae for non-Newtonian, power law fluids. The maximum shear rates 
(corrected for entrance effects) obtained in the rheometer with the current setup were in the 150,000 to 
170,000sec'' range. 

GFSSP simulations were performed with a flow circuit simulating the capillary rheometer and using 
Power Law gel viscosity coefficients from the experimental data. The agreement between the 
experimental data and the simulated flow curves was within ± 4% given quality entrance effect data. 
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Introduction 


The goal was to develop the experimental and 
analytical/computational tools required to predict 
the flow of non-Newtonian fluids through the 
various system components of a propulsion 
system: pipes, valves, pumps etc. To achieve this 
objective, we selected to augment the capabilities 
of NASA’s Generalized Fluid System Simulation 
Program (GFSSP) software. GFSSP is a general- 
purpose computer program designed to calculate 
steady-state and unsteady-state pressure and flow 
distributions in a complex fluid network. While 
the current version of the GFSSP code is able to 
handle various systems components, the implicit 
assumption in the code is that the fluid in the 
system is Newtonian. To extend the capability of 
the code to non-Newtonian fluids, such as silica 
gelled fuels and oxidizers, modifications to the 
momentum equations of the code were performed. 

We begin by examining the discretized momentum 
conservation equation shown below in Eq. 1. 
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Eq. 1 

This equation represents the balance of fluid 
forces acting on a given branch [Van Hooser, 
1999]. A typical branch configuration is shown 
below in Figure 1. As outlined in Van Hooser 
[1999], the two terms on the left side of Eq. 1 
represent the inertia of the fluid. The first term 
denotes the time dependence and must be 
considered for unsteady calculations. The second 
term is significant when there is a large change in 
area or density from branch to branch. The first 
term on the right side of the momentum equation 
represents the pressure gradient in the branch. The 
second term represents the effect of gravity. The 
gravity vector makes an angle (0) with the 
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assumed flow direction vector. The third term 
represents the frictional effect. Friction is 
modeled as a product of Kf, the square of the flow 
rate, and area. K f is a function of the fluid density 
in the branch and the nature of flow passage being 
modeled by the branch. The fourth term in the 
momentum equation represents the effect of the 
centrifugal force. This term will be present only 
when the branch is rotating as shown in Figure 1. 
K rot in this term is a factor representing fluid 
rotation and is unity when the fluid and the 
surrounding solid surface rotate at the same speed. 
This term also requires the radial distances from 
the upstream and downstream faces of the branch 
to the axis of rotation. 



Figure 1 - Schematic of a GFSSP Flow Branch 
Showing Sources of Momentum (from Van Hooser, 
1999) 

To account for non-Newtonian behavior, such as 
power law, the third term (friction effect) on the 
right side of the equation must be modified. The 
modification to the code is performed in a user 
subroutine, SORCEF, which is compiled and 
added to the GFSSP executable core at runtime. 
Figure 2 illustrates the interaction of the SORCEF 
user subroutine with the GFSSP solver module 
[Majumdar, 1999]. 
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Solver ModuleJWtion User Subrouti neW*™ 



Figure 2 - Interaction of User Subroutine with Solver Module (from Majumdar, 1999) 


The modifications to the friction momentum term 
(TERM4) in SORCEF consist of the calculation of 
a new shear rate for power law fluids and the 
relative friction force. 
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where n is the flow behavioral index, K is the 
consistency of the fluid, and r is the radius of the 
cross section of the flow branch under 
consideration. The implementation of the power 
law fluid behavior into the GFSSP code depends 
on knowledge of the two fluid coefficients, n and 
K. The determination of these parameters for the 
silica gels used in this program must be performed 
experimentally. The n and K parameters for silica 
water gels were determined experimentally at 
CFDRC’s Special Projects Laboratory using a 
capillary rheometer. The controller on the 
rheometer can be programmed to run in either 
controlled shear rate (constant flow-rate) or 
controlled shear stress (constant pressure) modes. 
The tests run with the silica gels were performed 


in constant shear rate mode. Data acquisition was 
performed with Lab VIEW. 

Silica Gel Testing 

The batch of silica gel used in all the tests was 
mixed with the CFDRC 10-gallon gel mixer. Prior 
to testing the gel was allowed to rest in the 
rheometer tank for at least twelve hours to ensure 
that the delicate structure of the gel had sufficient 
time to reform. During the tests silica gel was 
pressure fed and discharged through stainless steel 
pipes of l-36in length and three diameters 0.023, 
0.032, and 0.047in. The capillary tubes were 
soldered to l/4in nipples to be used with standard 
high-pressure fittings, and the tubes protruded out 
of the nipple 3/8 to l/2in. This resulted into a 
reentrant flow condition at the entrance of the 
capillary and an entrance pressure loss that needed 
to be evaluated experimentally. The pressure 
transducers were calibrated to within ±1% of 
reading by means of an Omega dead-weight tester. 
The flowrate was measured with an in-line 
McMillan turbine flowmeter. The accuracy of the 
flowmeter was checked by recording the time 
needed to fill calibrated laboratory beakers, and 
was found to be better than ± 2%. The 
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Figure 3 - Pressure Drop Versus Flowrate for 0.032in Tube of l-36in Length. Pressure Drop 

Includes Entrance Pressure Loss 


temperature of the gel during testing was 
monitored with an Omega DP-80 digital 
thermometer and was 21°C ± 1°C. The data 
collected in these tests were: pressure at tube 
entrance and volumetric flowrate. From these the 
uncorrected shear rate, shear stress, residence 
time, and viscosity were evaluated using formulae 
for non-Newtonian power law fluids. 


Tube 

Diameter 

rinl 

Length 

[in] 

Day # 

0.023 

34 

1 

0.032 

36 

1 

0.032 

36 

4 

0.047 

36 

1 (Test 1) 

0.047 

36 

1 (Test 2) 

0.047 

36 

4 (Test 1) 

0.047 

36 

4 (Test 2) 

0.032 

30 

5 

0.032 

24 

5 

0.032 

18 

5 

0.032 

12 

5 

0.032 

6 

5 

0.032 

1 

5 


Table 1 - Summary of Tests for Gelled Silica 


The maximum shear rates (corrected for entrance 
effects) obtained in the rheometer with the current 
setup were in the 150, 000-170, OOOsec' 1 . 

Table 1 lists the tests performed, the diameter 
/length of the stainless steel capillary tubes used, 
and the time passed between the gel mixing and 
the test time (in days). The tests were run over a 
one-week period in order to evaluate what the 
effect of gel settle-time is on the rheological 
properties of the gel. The 0.032in tube was cut in 
6in decrements and tested successively in day five 
in order to infer the entrance effects on gel flow. 

Analysis of the Gel Data 

In this section we present a preliminary analysis 
that has been performed on the gel data collected. 
The 0.032in tube data yielded the entrance loss 
coefficients for the reentrant flow. Figure 4.27 is a 
plot of the data collected in the six tests of day five 
and one test (36in) of day four. 

The entrance effect pressure drop was estimated 
according to the Couette method, as outlined in 
Ferguson and Kemblowski (1991). The pressure 
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drop per unit length of pipe according to this 
method is 

Eq . 3 

At i 2 -t, 

where AP 2 AP|, are measured pressure drops at the 
same flowrate, and L 2 , Li are the lengths of the 
tube tested. The method works best if the tubes 
are much longer than the entrance effect length. 
Thus, considering the data from the 30 in and 
24 in tests, the pressure drop per unit length of 
pipe was AP/AL = 38 psi/in. The pressure drop 
per unit length of pipe is varying slightly across 
the measured flowrate and tube length range, but 
for simplicity the 38 psi/in value will be used to 
estimate the pressure loss due to reentrant flow at 
tube entrance. 

Choosing from the middle of the flowrate/length 
range, 600 mL/min and 18 in, the estimated 
pressure loss due to entrance effects is, 

AP e = APtot - (AP/AL)-L = 824psi - 38 psi/in- 18 in 
= 140 psi 

It is useful to define the pressure loss effect in 
terms of the dynamic pressure in the tube 
[Munson, 1990]. 


K 


e 





Eq. 4 


For the 18 in tube flowing 600 mL/min the 
average velocity is V = 19.4 m/sec, and the 
dynamic pressure is 27.12 psi. The entrance loss 
coefficient is K e = 5.16. We note that the pressure 
loss due to entrance effect is five times greater 
than the dynamic pressure in the tube. In 
comparison, the reentrant loss coefficient for water 
flow [Munson, 1990] is only K e v/a ter = 0.8. 
Generally, the entrance loss coefficient is 
dependent on geometry and Reynolds number. 
However, in most cases of interest the loss 
coefficient is only a function of geometry. For 


now we will assume a constant entrance loss 
coefficient for all tubes of K e = 5 (It will be shown 
in the next section that this assumption introduces 
some error in the determination of the flow-curves 
for different pipe diameters). Using this value and 
the average velocity, the corrected pressure drop 
along the tubes can be calculated. 

The non-Newtonian, power law fluid shear stress 
and shear rate are described by the following 
equations, 
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Eq. 5 


where L is the length of the tube, d is the diameter, 
k is a consistency factor (Pa-sec"), and n is the 
behavioral index. The flow curve of a power law, 
time independent fluid is 


v = k-y n 


Eq. 6 


The determination of k and n can be done by 
taking advantage of the quasi-linear behavior of 
the data in log-log representation. Replacing the 
shear stress and shear rate with the values in Eq. 5, 
the shear stress becomes. 
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Eq.7 


We solve for k and n by fitting Eq. 7 to the data 
and obtain n = 0.349 and k = 23.06 Pa-sec". The 
quasi-linear behavior of the data scaled 
logarithmically is evident in Figure 4. 
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Figure 4 - Logarithmic Representation of the tvs. y Data 


Thus Eq. 6 is rewritten as 

r = 23.06 -y 8349 Eq. 8 

These values are consistent with data obtained in 
previous experiments with a stainless steel tube of 
0.084in diameter and 6ft length. Additionally, the 
behavioral index of the silica gel is very similar to 
the behavior index of buttermilk, as reported by 
Butler and O’Donnell (1999). The viscosity of the 
silica gel is therefore 

p = 23.06 -y °- 349-1 Eq. 9 

The flow curve of the silica gel is also time 
dependent. Generally, this dependence is 
expressed by multiplying Eq. 9 by a time 
dependent function [Butler, 1999], 

^,r) = ^(0-[ 23 - 06 -r°' 349 “ 1 ] Eq. 10 

where, in general X(t) is a function of time. As 
specified by Butler, (1999) and Eq. 10, the 
viscosity of the time dependent shear thinning gel 
can be expressed as a product between a unitless. 


time dependent function X(t), and a power law, as 
determined in Eq. 9. The time dependent function, 
7,(t) tracks in fact the evolution of a structural 
parameter X, which is described by a second order 
kinetic equation derived by Petrellis and 
Flumerfelt, (1973), 

^■ = -k l -(X-A e f, Eq. 11 


for X >X e . The structural parameter X ranges from 
an initial value of unity (at zero time of shear) to 
an equilibrium value (A, e <1 ) for t -»oo and k t is a 
rate-constant, which is a function of shear rate. 
Integrating Eq. 1 1 we obtain, 


d[X-X ) 

— (.i — A )~ = —k. • t — C , , 


X(f) = X 


e + C 1 +k l -t 


Eq. 12 
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where Ci is a constant of integration. The 
determination of A, e , ki, and Q was done by 
utilizing the already available shear thinning (i.e. 
power law behavior) data and by employing our 
Brookfield DV-II+ programmable viscometer with 
WinGather™ software to determine the time 
variation of silica gel viscosity at constant shear 
rate. A set of five LV-series Brookfield spindles 
was used to make viscosity measurements over a 
five order-of-magnitude (1 O' 4 - 10 Pa- sec) range. 

Extrapolating the measured data, the zero shear 
viscosity for the water silica gel tested (12.5% 
silica by mass) was 8 Pa-sec. At the limited shear 
rate of 7-10 s' 1 exerted by the spindle of the 
viscometer, after 1.5 hours, the viscosity of the gel 
settled at around 0.2 Pa- sec and did not change 
further to a noticeable degree. Moreover, unlike 
the gel sheared by the high shear rheometer, this 
gel recovered its original viscosity once allowed to 
sit in the test cup for about half an hour. 

The estimation of the three unknown constants in 
Eq. 12 begins by noticing that at the limit of t -»oo 
the structural parameter Aft) — > A, e . At t = 4000 
sec, the viscosity of the gel was 0.2 Pa-sec and 
virtually constant. Thus we can write 

f e -k-y"- x =0.2Pa-SQC Eq. 13 

with k = 23.06 Pa-sec", dy/dt = 7.3 sec" 1 , and n = 
0.349. Solving, we obtain A. e = 0.03. Also, we 
note that at t = 0, the viscosity is 8 Pa-sec and Eq. 
12 reduces to, 
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Eq. 14 


Eq. 15 


Solving for Q we obtain Q = 0.81. To obtain ki 
we simply remove the shear thinning contribution 
from the data set. These data is then fitted with a 
function of the form, 


A + 0.03 = 


1 

0.81 + k\ • t 


Eq. 16 


from where the value of ki = 0.022sec"’ is 
obtained. Summarizing, the silica gel viscosity 
can thus be written as. 


f 


p{t,y)= 0.03 + — , 

l 0.81 + 0.022- sec - • 

• (23.06.Pa • sec" )• {f) nX . 


Eq. 17 


Eq. 17 fits the original data within ±25% as shown 
in Figure 5. 



Time (sec) 


Figure 5 - Original Viscosity Data and Eq. 1 7 Fit 
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P vs. Q Silica Water Gel Flow Data and GFSSP Simulation 
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Figure 6 - Experimental Silica Gel Data and GFSSP Simulations 


Modeling Results 

The GFSSP simulation included a simple flow 
circuit similar to the capillary rheometer (see 
Figure 7). 



Figure 7 - GFSSP Schematic of Capillary 
Rheometer 


The boundary nodes specify the tank and 
discharge pressure and the parameters of the pipe 
section are set to match the actual pipes tested. 
The test data for the 0.032 in and 0.047 in 
diameter tubes, and the corresponding simulated 
flow curves are shown in Figure 6. 


We note from Figure 6 that the 0.032 in tube data 
is approximated much better (±4%) than the 
0.047 in tube data (±15%). Since the n and K 
values were calculated with the 0.047 in tube data, 
and the entrance loss coefficient was calculated 
from the 0.032 in data, we can conclude that the 
error in the 0.047 in simulation is due to an 
overestimation of the pressure loss at the entrance 
of the 0.047 in tube. Also, we can conclude that 
the relative insensitivity of the n and K 
determination to the error made in the evaluation 
of the entrance effect is due to the fact that the 
entrance loss manifests itself as a uniform bias of 
the pressure data, and does not greatly affect the 
shape of the functional relationship expressed by 
Eq. 6. These considerations will permit the 
estimation of the entrance effect from the GFSSP, 
if some experimental pressure/flowrate data is 
available. 
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Conclusions 

In this study we successfully implemented into 
GFSSP flow equations for fluids with power law 
behavior. This required the modification of the 
sink terms of the momentum equation. This was 
the first step in the introduction of non-Newtonian, 
time dependent flow treatment into the GFSSP 
code. As part of the study, we also developed and 
experimentally verified the procedure for 
quantifying the properties of non-Newtonian, 
time-dependent fluids necessary to define the 
physical properties of the fluids. Future efforts 
will be directed towards expanding the code to 
incorporate treatment of the time dependent 
behavior. 
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